#include "analysis.h"
#include "utilities.h"
#include "rngs.h"
#include<stdio.h>
#include<stdlib.h>
#include<unistd.h>
#include<sys/wait.h>
#include<string.h>

/*--------------------------------------------*/
void init_guess(int num_ensemble)
{
	int status1 = 0;
	int info = 0;
	int i = 0;
	int j = 0;
	char strfile[200];
	char command[300];

	// Operate on server to get configuration ready
	if (fork())
	{
		wait(&status1);
	}
	else
	{
		info = system("ssh ITE366CS@130.85.94.246 'cd /cygdrive/c/hysplit4/working/; cp exp3/* .'");
		exit(0);
	}			
	// Copy CONTROL file
	if (fork())
	{
		wait(&status1);
	}
	else
	{
		info = system("scp ITE366CS@130.85.94.246:/cygdrive/c/hysplit4/working/CONTROL /home/shiming1/InfoSphereStreams/letkf-spade/hysplit/data > /dev/null");
		exit(0);
	}

	modify_control();
	// Copy k
	if (fork())
		{
		wait(&status1);
		}
		else
		{
		info = system("ssh ITE366CS@130.85.94.246 'cd /cygdrive/c/hysplit4/working/; ./par2asc.exe -iPARDUMP -otes; ../exec/./con2asc -s -icdump -otestn; mv PARDUMP PARINIT; ../exec/./hycs_std; cp PARDUMP PARDUMPk; cp cdump cdumpk'");
		exit(0);
		}	
	
	//Generate num_ensemble initial guess, running Hysplit multi-times
	//1. copy tes
	if (fork())
	{
		wait(&status1);
	}
	else
	{
		info = system("scp ITE366CS@130.85.94.246:/cygdrive/c/hysplit4/working/tes /home/shiming1/InfoSphereStreams/letkf-spade/hysplit/data > /dev/null");
		exit(0);
	}
	//2.Read into array and make purturbations
	int numpar = 0;
	char str1[100];
	char str2[100];
	char str3[100];
	char str4[100];
	char str5[100];
	double d_rand1 = 0.0;
	double d_rand2 = 0.0;
	FILE *fp;
	fp = fopen("/home/shiming1/InfoSphereStreams/letkf-spade/hysplit/data/tes","r");
	fgets(str1, 100, fp);
	fgets(str2, 100, fp);
	fgets(str3, 100, fp);
	fgets(str4, 100, fp);
	info = fscanf(fp,"%d", &numpar);
	fgets(str5, 100, fp);
	double *mass = (double *)malloc(numpar*sizeof(double) );
	float *tlat = (float *)malloc(numpar*sizeof(float) );
	float *tlon = (float *)malloc(numpar*sizeof(float) );
	float *zlvl = (float *)malloc(numpar*sizeof(float) );
	float *sigh = (float *)malloc(numpar*sizeof(float) );
	float *sigw = (float *)malloc(numpar*sizeof(float) );
	float *sigv = (float *)malloc(numpar*sizeof(float) );
	int16_t *page = (int16_t *)malloc(numpar*sizeof(int16_t) );
	int16_t *hdwp = (int16_t *)malloc(numpar*sizeof(int16_t) );
	int16_t *ptyp = (int16_t *)malloc(numpar*sizeof(int16_t) );
	int16_t *pgrd = (int16_t *)malloc(numpar*sizeof(int16_t) );

	for(i=0; i<numpar; i++)
	{
		info = fscanf(fp,"%lE", mass+i);
		info = fscanf(fp,"%f", tlat+i);
		info = fscanf(fp,"%f", tlon+i);
		info = fscanf(fp,"%f", zlvl+i);
		info = fscanf(fp,"%f", sigh+i);
		info = fscanf(fp,"%f", sigw+i);
		info = fscanf(fp,"%f", sigv+i);
		info = fscanf(fp,"%hd", page+i);
		info = fscanf(fp,"%hd", hdwp+i);
		info = fscanf(fp,"%hd", ptyp+i);
		info = fscanf(fp,"%hd", pgrd+i);
		info = fscanf(fp,"%d", &info);
	}
	fclose(fp);

	for(i=0; i<num_ensemble; i++)
	{
		fp = fopen("/home/shiming1/InfoSphereStreams/letkf-spade/hysplit/data/parinittxt","w");

		fprintf(fp,"%s", str1);
		fprintf(fp,"%s", str2);
		fprintf(fp,"%s", str3);
		fprintf(fp,"%s", str4);
		fprintf(fp,"%12d", numpar);
		fprintf(fp,"%s", str5);
	
		for(j=0; j<numpar; j++)
		{
			d_rand1 = Random();
			d_rand2 = Random();
			fprintf(fp,"%12.3E\n", mass[j]+2.0*(d_rand1-d_rand2) );
			fprintf(fp,"%12.3f", tlat[j]);
			fprintf(fp,"%12.3f", tlon[j]);
			fprintf(fp,"%12.3f", zlvl[j]);
			fprintf(fp,"%12.3f", sigh[j]);
			fprintf(fp,"%12.3f", sigw[j]);
			fprintf(fp,"%12.3f\n", sigv[j]);
			fprintf(fp,"%12d", page[j]);
			fprintf(fp,"%12d", hdwp[j]);
			fprintf(fp,"%12d", ptyp[j]);
			fprintf(fp,"%12d", pgrd[j]);
			fprintf(fp,"%12d\n", j+1);
		}
		fclose(fp);
	// Copy back and convert to PARINIT and run for cdump
	if (fork())
	{
		wait(&status1);
	}
	else
	{
		info = system("scp /home/shiming1/InfoSphereStreams/letkf-spade/hysplit/data/parinittxt ITE366CS@130.85.94.246:/cygdrive/c/hysplit4/working/ > /dev/null");
		exit(0);
	}

	if (fork())
	{
		wait(&status1);
	}
	else
	{
		sprintf(command,"ssh ITE366CS@130.85.94.246 'cd /cygdrive/c/hysplit4/working/; ./asc2par.exe -iparinittxt -oPARINIT; ../exec/./hycs_std; ./par2asc.exe -iPARDUMP -otes; mv tes tes%d; ../exec/./con2asc -s -icdump -otest; mv test.txt test%d; '", i, i);
		info = system(command);
		exit(0);
	}
	// convert cdump tp test%d copy back
	if (fork())
	{
		wait(&status1);
	}
	else
	{
		sprintf(command,"scp ITE366CS@130.85.94.246:/cygdrive/c/hysplit4/working/tes%d /home/shiming1/InfoSphereStreams/letkf-spade/hysplit/data", i);
		info = system(command);
		exit(0);
	}

	if(fork())
	{
		wait(&status1);
	}
	else
	{
		sprintf(command,"scp ITE366CS@130.85.94.246:/cygdrive/c/hysplit4/working/test%d /home/shiming1/InfoSphereStreams/letkf-spade/hysplit/data", i);
		info = system(command);
		exit(0);
	}

  }// end of for_loop i 

}
/*-------------------------------------------------------------------------*/
int single_time_step(int nstep, int astep, int num_ensemble)

{
// Write xb_bar_g to PARINIT
	int status1;
	int info = 0;
	int i = 0;
	char command0[150]="ssh ITE366CS@130.85.94.246 'cd /cygdrive/c/hysplit4/working/; ./asc2par.exe -i";
	char cmd0bak[150]="ssh ITE366CS@130.85.94.246 'cd /cygdrive/c/hysplit4/working/; ./asc2par.exe -i";
	char command1[150];
	char command2[300]="ssh ITE366CS@130.85.94.246 'cd /cygdrive/c/hysplit4/working/;";
	char cmd2bak[300]="ssh ITE366CS@130.85.94.246 'cd /cygdrive/c/hysplit4/working/;";
	char pardump[100];
	char cdump[100];
	char parinit[20];
	
	if ( nstep%astep == 0)
		{
			// Copy it back to server
			if (fork())
			{
			wait(&status1);
			}
			else
			{
			info = system("ssh ITE366CS@130.85.94.246 'cd /cygdrive/c/hysplit4/working/; ./asc2par.exe -iparinittxt -oPARINIT'");
			exit(0);
			}
		}
	else
		{
		if (fork())
			{
			wait(&status1);
			}
		else
			{
			info = system("ssh ITE366CS@130.85.94.246 'cd /cygdrive/c/hysplit4/working/; cp PARDUMPk PARINIT'");
			exit(0);
			}	
		}			
		// Change the CONF and CONTROL files
		modify_control();
	
		// Run the hycs_std.exe

		if (fork())
		{
		wait(&status1);
		}
		else
		{
		info = system("ssh ITE366CS@130.85.94.246 'cd /cygdrive/c/hysplit4/working/; ../exec/./hycs_std'");
		exit(0);
		}
		char command3[300];
		sprintf(command3, "ssh ITE366CS@130.85.94.246 'cd /cygdrive/c/hysplit4/working/; mv PARDUMP PARDUMPk; cp cdump bakdump%d; mv cdump cdumpk'", nstep);
		if (fork())
		{
		wait(&status1);
		}
		else
		{
		info =system(command3);

		exit(0);
		}
	// Also for a cdump copy of 10000m height, which just rerun it.

	// About the H_en
	// Write each line to one PARINIT file
	for (i=0; i<num_ensemble; i++)
	{
		// Convert parinittxtn to PARINIT ./asc2par -iparinittxt -oPARINIT

		sprintf(parinit,"parinittxt%d",i);
		strcpy(command0,cmd0bak);
		strcat(command0, parinit);
		strcat(command0, " -oPARINIT'"); 
		// Copy PARINIT_n to PARINIT, copy CONTROLBAK CONTROL
		if (fork())
		{
		wait(&status1);
		}
		else
		{
			if ( ((nstep)%(astep) == 0) )//|| (nstep == 1))
				{
				info = system(command0);
				}
			else
				{
				sprintf(pardump,"cp PARDUMP%d PARINIT; '",i);
				strcpy(command1, cmd2bak);
				strcat(command1, pardump);
				info = system(command1);
				}	
		exit(0);
		}
		
		// Run hycs_std
		if (fork())
		{
		wait(&status1);
		}
		else
		{
		info = system("ssh ITE366CS@130.85.94.246 'cd /cygdrive/c/hysplit4/working/; ../exec/./hycs_std'");// cp CONTROLBAK CONTROL; 
		exit(0);
		}

		// Copy PARDUMP to PARDUMP_n, copy cdump to cdump_n
		
		if (fork())
		{
		wait(&status1);
		}
		else
		{
		sprintf(pardump,"./par2asc.exe -iPARDUMP -otes%d; mv PARDUMP PARDUMP%d; ", i, i );
		sprintf(cdump,"../exec/./con2asc -s -icdump -otestn; mv testn.txt test%d; mv cdump cdump%d;'", i, i );
		strcpy(command2, cmd2bak);
		strcat(command2,pardump);strcat(command2,cdump);
		info = system(command2);
		exit(0);
		}

		// Copy tes%d test%d to calculation side
		if (fork())
		{
			wait(&status1);
		}
		else
		{
		sprintf(command3, "scp ITE366CS@130.85.94.246:/cygdrive/c/hysplit4/working/tes%d /home/shiming1/InfoSphereStreams/letkf-spade/hysplit/data/", i );
		info = system(command3);
		exit(0);
		}

		if (fork())
		{
		wait(&status1);
		}
		else
		{
		sprintf(command3, "scp ITE366CS@130.85.94.246:/cygdrive/c/hysplit4/working/test%d /home/shiming1/InfoSphereStreams/letkf-spade/hysplit/data/", i );
		info = system(command3);
		exit(0);
		}
	
	} // end of for_loop num_ensemble

return 0;
}

/*-----------------------------------------------------------------------------*/
void all_dist(float *dist_pardump, float *dist_cdump, float *dist_datem, 
	float *xb_tlat, float *xb_tlon, float *xb_zlvl, float *yb_clat, float *yb_clon, float *obs_lat, float *obs_lon, int numpar, int numline, int numobs)
{
int i = 0;
int j = 0;
int shift = 0;

for (i=0; i<numpar; i++)
	{
	for (j=(i+1); j<numpar; j++)
		dist_pardump[shift+j-1-i] = (xb_tlat[i]-xb_tlat[j])
									*(xb_tlat[i]-xb_tlat[j])
									+(xb_tlon[i]-xb_tlon[j])
									*(xb_tlon[i]-xb_tlon[j])
									+(xb_zlvl[i]-xb_zlvl[j])
									*(xb_zlvl[i]-xb_zlvl[j]);

	shift = shift + numpar - i - 1;
	}	
for (i=0; i<numline; i++)
	for (j=0; j<numpar; j++)
		dist_cdump[i*numpar+j] = (xb_tlat[j]-yb_clat[i])
								*(xb_tlat[j]-yb_clat[i])
								+(xb_tlon[j]-yb_clon[i])
								*(xb_tlon[j]-yb_clon[i]);

for (i=0; i<numpar; i++)
	for (j=0; j<numobs; j++)
		dist_datem[i*numobs+j] = (xb_tlat[i]-obs_lat[j])
								*(xb_tlat[i]-obs_lat[j])
								+(xb_tlon[i]-obs_lon[j])
								*(xb_tlon[i]-obs_lon[j]);		

}

/*-----------------------------------------------------------------------*/
int find_neighbors(int i, int radius, int num_ensemble, int numpar, int numline, int numobs, double *xb_bar, double *yb_bar, double **en_mass, double **en_con, double *obs, double *localxb_bar, double *localyb_bar, double **localXb, double **localYb, double *localobs, double *localR, int *len_r, float *dist_pardump, float *dist_cdump, float *dist_datem, double valid_dist)
{
// Find the smallest radius number of elements, and store in l_state_p
int k = 0;
int j = 0;
int shift = 0;
int l_len = radius;
int xpar = 0;
int *xpos = (int *)malloc(radius*sizeof(int));
int *ypos = (int *)malloc(radius*sizeof(int));
float *l_dist = (float *)malloc(radius*sizeof(float));
for (k=0; k<radius; k++)
	{
	l_dist[k] = 0.0;
	xpos[k] = 0.0;
	ypos[k] = 0.0;
	}

//get 100m level points
shift = i * numobs;
for (j=0; j<radius; j++)
	{
		l_dist[j] = dist_datem[shift+j];
		ypos[j] = j;
	}
small_sort_dist(l_dist, ypos, radius);
for (j=radius; j<numobs; j++)
	{
		for (k=0; k<radius; k++)
		{
		if (dist_datem[shift+j] < l_dist[k])
			{
				l_dist[k] = dist_datem[shift+j];
				ypos[k] = j;
				break;
			}
		}
	}
//--------------------
for (j=0; j<radius; j++)
	{
	if (l_dist[j] > valid_dist)
		l_len--;
	}
//--------------------
for (j=0; j<l_len; j++)
	{
	localobs[j] = obs[ypos[j]];
	}

// Collect Xb_g, xb_abr local
for (j=0; j<l_len; j++)
	{
	xpar = find_min_dist(dist_datem, ypos[j], numpar, numobs);
	localxb_bar[j] = xb_bar[xpar];
	xpos[j] = xpar;
	for (k=0; k<num_ensemble; k++)
		{
		localXb[k][j] = en_mass[k][xpar];
		}
	}
// Collect Yb_g, yb_bar
for (j=0; j<l_len; j++)
	{
	xpar = find_min_dist(dist_cdump, xpos[j], numline, numpar);
	localyb_bar[j] = yb_bar[xpar];
	for (k=0; k<num_ensemble; k++)
		{
		localYb[k][j] = en_con[k][xpar];
		}
	}

// l_state_p->R
if (l_len > 0)
	{
		localR = (double *)realloc(localR, l_len*l_len*sizeof(double));
		for(j=0; j<l_len*l_len; j++ )
			localR[j] = 0.0;

		for (k=0; k<l_len; k++)
			{
			localR[k+k*l_len] = 0.001;
			}
	}
*len_r = l_len;

free(xpos); 
free(ypos);
free(l_dist);

return 0;
}

/*-----------------------------------------------------------------*/
int find_min_dist(float *dist_array, int pos, int sizex, int sizey)
{
int i = 0;
int resx = 0;
float tmp = 999999.0;
for (i=0; i< sizex; i++)
	{
	if (tmp > dist_array[i*sizey+pos])
		{
		tmp = dist_array[i*sizey+pos];
		resx = i;
		}
	}
	return resx;
}

/*------------------------------------------------------------*/
void small_sort_dist(float *numbers, int *pos, int array_size)
{// heap sort
  int i, tempos;
  float temp;
 
  for (i = array_size - 1; i >= 0; i--)
    siftDown(numbers, pos, i, array_size - 1);
 
  for (i = array_size-1; i >= 1; i--)
  {
    temp = numbers[0];
    numbers[0] = numbers[i];
    numbers[i] = temp;
    tempos = pos[0];
    pos[0] = pos[i];
    pos[i] = tempos;
    siftDown(numbers, pos, 0, i-1);
  }
}

/*----------------------------------------------------------*/
void siftDown(float *numbers, int *pos, int root, int bottom)
{
  int done, maxChild, tempos;
  float temp;
  
  done = 0;
  while ((root*2 <= bottom) && (!done))
  {
    if (root*2 == bottom)
      maxChild = root * 2;
    else if (numbers[root * 2] > numbers[root * 2 + 1])
      maxChild = root * 2;
    else
      maxChild = root * 2 + 1;
 
    if (numbers[root] < numbers[maxChild])
    {
      temp = numbers[root];
      numbers[root] = numbers[maxChild];
      numbers[maxChild] = temp;
      root = maxChild;
      tempos = pos[root];
      pos[root] = pos[maxChild];
      pos[maxChild] = tempos;
    }
    else
      done = 1;
  }
}
